#pragma once

//Rotate left and right by (1+-ig5)/sqrt(2)
//We distinguish for types of rotations, ir=rsink*2+rsource
//We divide the spinspin into 4 blocks, according to id<2
// -------------------------------
// | div |  0  |  1  |  2  |  3  |
// -------------------------------
// | 0 1 | + 1 | 1 + | 1 - | - 1 |
// | 2 3 | 1 - | - 1 | + 1 | 1 + |
// -------------------------------
// so just have to specify which is the block which rotate as +-i
// From Francesco
void rotate_spinspin_to_physical_basis(spinspin s,int rsi,int rso)
{
  const int list_prb[4]={0,1,2,3},list_mrb[4]={3,2,1,0}; //plus and minus rotating blocks
  const int so_shft[4]={0,2,0,2},si_shft[4]={0,0,2,2};   //start of dirac indexes defining blocks

  int ir=rsi*2+rso,prb=list_prb[ir],mrb=list_mrb[ir];

  for(int dso=0;dso<2;dso++)
    for(int dsi=0;dsi<2;dsi++)
      {
        int pso=dso+so_shft[prb],psi=dsi+si_shft[prb];
        int mso=dso+so_shft[mrb],msi=dsi+si_shft[mrb];

        // rotation with +,-
        s[pso][psi]=I*s[pso][psi];
        s[mso][msi]=-I*s[mso][msi];
      }
}

void rotate_colorcolorspinspin_to_physical_basis(colorcolorspinspin S,int r)
{
  printf("Rotating ...\n");

  for (int ic1=0; ic1<3; ic1++)
     for (int ic2=0; ic2<3; ic2++) 
     {
              spinspin Saux;
              for (int is1=0; is1<4; is1++)
                  for (int is2=0; is2<4; is2++)
                           Saux[is2][is1]=S[ic2][ic1][is2][is1];
              rotate_spinspin_to_physical_basis(Saux,r,r);
              for (int is1=0; is1<4; is1++)
                    for (int is2=0; is2<4; is2++)
                            S[ic2][ic1][is2][is1]=Saux[is2][is1];
    }

}

/*
void complex_prod(complex a,complex b,complex c)
{
	a=b*c;
}


void rotate_spinspin_to_physical_basis(spinspin out, int r)
{

  if (r==1){
  //Pplus*s
  spinspin temp;
  memset(temp,0,sizeof(spinspin));

  for(int isfila=0;isfila<4;isfila++)
    for(int iscolumn=0;iscolumn<4;iscolumn++)
    {
        complex_prod(temp[isfila][iscolumn],Pplus.entr[iscolumn],out[isfila][Pplus.pos[iscolumn]]);
    }
  memcpy(out,temp,sizeof(spinspin));

  //(Pplus*s)*Pplus
  memset(temp,0,sizeof(spinspin));
  for(int isfila=0;isfila<4;isfila++)
    for(int iscolumn=0;iscolumn<4;iscolumn++)
    {
        complex_prod(temp[isfila][iscolumn],out[isfila][Pplus.pos[iscolumn]],Pplus.entr[iscolumn]);
    }
  memcpy(out,temp,sizeof(spinspin));
 }

  if (r==0){
  //Pminus*s
  spinspin temp;
  memset(temp,0,sizeof(spinspin));

  for(int isfila=0;isfila<4;isfila++)
    for(int iscolumn=0;iscolumn<4;iscolumn++)
    {
        complex_prod(temp[isfila][iscolumn],Pminus.entr[iscolumn],out[isfila][Pminus.pos[iscolumn]]);
    }
  memcpy(out,temp,sizeof(spinspin));

  //(Pminus*s)*Pminus
  memset(temp,0,sizeof(spinspin));
  for(int isfila=0;isfila<4;isfila++)
    for(int iscolumn=0;iscolumn<4;iscolumn++)
    {
        complex_prod(temp[isfila][iscolumn],out[isfila][Pminus.pos[iscolumn]],Pminus.entr[iscolumn]);
    }
  memcpy(out,temp,sizeof(spinspin));

  }
}
*/
